%针对N>13的原子体系，在球坐标系中进行构型空间采样
%初始化参数
N = 55;
d = 0.88; 

%导入N=13的构型
DATA = load('.\SA_revised\N=13.txt');
data = DATA(15:27,:);
%en_13 = -11.08158890;

%初始化构型
config0 = cell(1,N);
for i = 1:13
    config0{i} = data(i,:);
end

%在N=13的基础上在外层随机产生一个构型，config_sph_former从N=14开始是用球坐标表示
[config_sph_former, en_former] = create_config(config0);

%en_min记录退火过程中目前最小能量对应的构型{1}和能量{2}
en_min = cell(1,2);
en_min{1} = config_sph_former;
en_min{2} = en_former;

T_max = 5;
T_min = 0.5;
lambda = 0.92;
step = 20000;
dr = 0.5;
dtheta = pi;
dphi = 2*pi;
T = T_max;

while T > T_min

    for i = 1:step
        [config_sph,en] = change_config_sph(config_sph_former, dr, dtheta, dphi);
        if(en < en_former)
            config_sph_former = config_sph;
            en_former = en;
            if(en < en_min{2})
                en_min{1} = config_sph;
                en_min{2} = en;
            end
        else
            rand_num = rand;
            factor = exp((en_former-en)/T);
            if(rand_num < factor)
                config_sph_former = config_sph;
                en_former = en;
            else

            end
        end
    end
    T = lambda * T;

end

%将球坐标转化为直角坐标
config_min = sph2car(en_min{1});
en_min{1} = config_min;
config_former = config_min;
en_former = en_min{2};

T_max = 5 * 1e-1;
T_min = 1 * 1e-3;
lambda = 0.9;
step = N*100;
step_r = 0.1;
choice = 0;
limit = [-3,3,-3,3,-3,3];
T = T_max;

while T > T_min
    for i = 1:step
        [config,en] = change_config(config_former, en_former, limit,step_r,choice);
        if(en < en_former)
            config_former = config;
            en_former = en;
            if(en < en_min{2})
                en_min{1} = config;
                en_min{2} = en;
            end
        else
            rand_num = rand;
            factor = exp((en_former-en)/T);
            if(rand_num < factor)
                config_former = config;
                en_former = en;
            else

            end
        end
    end
    T = lambda * T;
    step_r = lambda^(5/7) * step_r;

end

%写.txt文件，存数据
fname = sprintf('N=%d.txt',N);
fid = fopen(['.\SA_revised\',fname],'w');
for k = 1:N
    for m = 1:3
        fprintf(fid, '%.8f\t', config_former{k}(m));
    end
    fprintf(fid, '\n');
end
for m = 1:3
    fprintf(fid, '%.8f\t', en_former);
end
fprintf(fid, '\n');
for k = 1:N
    for m = 1:3
        fprintf(fid, '%.8f\t', en_min{1}{k}(m));
    end
    fprintf(fid, '\n');
end
for m = 1:3
    fprintf(fid, '%.8f\t', en_min{2});
end
fclose(fid);

